# ========== 第一部分：方差分析 ==========
# 设置工作目录
setwd("D:/Rassignment/第三次作业")

# 加载包
library(openxlsx)
library(tidyverse)

# 读取数据
ADdata <- read.xlsx("merged_data.xlsx") %>%
  column_to_rownames("X")  # 设置行名为X列

# 数据处理
AD <- log2(ADdata) %>%
  replace(is.infinite(.), 0) %>%  # 将无穷值替换为0
  select(order(colnames(.)))  # 按列名排序

# 转置数据并添加分组信息
AD_long <- AD %>%
  t() %>%
  as.data.frame() %>%
  rownames_to_column("sample") %>%
  mutate(
    group = str_replace_all(sample, "\\d", ""),  # 简化样本名
    group = factor(group)  # 转换为因子
  ) %>%
  select(-sample)  # 移除样本名列

# 提取基因名
genenames <- colnames(AD_long)[-ncol(AD_long)]

# 创建空向量存储p值
p_values <- setNames(rep(NA, length(genenames)), genenames)

# 遍历基因进行方差分析
for (gene in genenames) {
  # 筛选非零值样本
  valid_samples <- AD_long %>%
    filter(!!sym(gene) != 0) %>%
    group_by(group) %>%
    summarise(n = n(), .groups = "drop")
  
  # 检查每组至少3个样本
  if (all(valid_samples$n >= 3)) {
    # 构建公式
    formula <- as.formula(paste(gene, "~ group"))
    
    # 执行方差分析
    anova_model <- aov(formula, data = AD_long)
    
    # 提取p值
    p_value <- summary(anova_model)[[1]]$`Pr(>F)`[1]
    p_values[gene] <- p_value
  }
}

# 转换为数据框
p_values_df <- data.frame(P = p_values, row.names = names(p_values))

# ========== 第二部分：GO富集分析 ==========
# 设置工作目录
setwd("D:/Rassignment/第四次作业")

# 加载包
library(clusterProfiler)
library(org.Hs.eg.db)
library(ggplot2)

# 加载数据
load("volcano.RData")  # 加载prostat对象

# 筛选显著蛋白
sigpro <- subset(prostat, P < 0.05)
gene_list <- sigpro$ID

# GO富集分析
go_result <- enrichGO(
  gene = gene_list,
  OrgDb = org.Hs.eg.db,
  keyType = "SYMBOL",
  ont = "BP",  # 生物过程
  pAdjustMethod = "BH",  # Benjamini-Hochberg校正
  pvalueCutoff = 0.05,
  qvalueCutoff = 0.05,
  readable = FALSE
)

# 结果处理
go_result_df <- as.data.frame(go_result) %>%
  arrange(p.adjust) %>%  # 按校正p值排序
  head(30)  # 取前30个显著结果